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Molecular dynamics (MD) simulation with modified Brenner's reactive empirical bond order (REBO) po- 
tential is a powerful tool to investigate plasma wall interaction on divertor plates in a nuclear fusion device. 
However, MD simulation box's size is less than several nm for the performance of a computer. To extend the 
size of the MD simulation, we develop a hybrid simulation code between MD code using REBO potential and 
binary collision approximation (BCA) code. Using the BCA code instead of computing all particles with a high 
kinetic energy for every step in the MD simulation, considerable computation time is saved. By demonstrating 
a hydrogen atom injection on a graphite by the hybrid simulation code, it is found that the hybrid simulation 
code works efficiently in a large simulation box. 



I. INTRODUCTION 

To understand a mechanism of chemical and physical in- 
teractions between hydrogen plasmas and divertor plates in a 
nuclear fusion device, it is necessary to study an elementary 
processes of the reactions. In molecular dynamics (MD) sim- 
ulations, the equations of motion of atoms are solved numeri- 
cally. We investigated plasma- wall interaction on the divertor 
plates made of carbon materials by the MD simulation with 
modified Brenner's reactive empirical bond order (REBO) po- 
tential^'^ in the previous works^"^. A typical scale length of 
the MD simulation box in these works are several nm. 

In order to expand the simulation box to more realistic scale 
length, i.e., several /im, we develop a hybrid simulation code 
between the MD simulation code with the REBO potential 
and atomic collision in any structured target(ACVT) code^. 
In the ACVT code, which uses binary collision approxima- 
tion (BCA), a two-body interaction is calculated instead of 
computing all particles for every step in the MD simulation. 
Therefore computation time is saved. We choose the ACVT 
code as the partner of the MD simulation for its fast calcu- 
lation of a location and a velocity of the particle. The main 
problem of hybridization between the two simulation codes is 
how to define a threshold kinetic energy of the injection par- 
ticle. Eor higher kinetic energy of particle than the threshold 
kinetic energy, we calculate the trajectory of the particle by 
the ACVT code. As the kinetic energy is dissipated to its sur- 
roundings and attaches to the threshold kinetic energy £^th, 
we give a position and a velocity of the projectile to the MD 
simulation code as the initial condition. Then we perform the 
MD simulation. Generally speaking, because a multi-body in- 
teraction with the neighborhood of the projectile must be cal- 
culated in the MD simulation code with the REBO potential, 
computation time becomes longer than in the ACVT code. To 
make the performance of the simulation better, it is necessary 
to shorten the MD simulation and extend the ACVT simulation 
box. However, it is known that the ACVT code cannot calcu- 
late exactly the particle motion in a low energy. Comparing 
the MD simulation with the ACVT simulation, two simula- 
tion results are consistent in higher kinetic energy than 200 



eV— . Thus we set £^th=200 eV at which a simulation algo- 
rithm changes from the ACVT simulation to the MD simula- 
tion. 



II. ACVT-MD HYBRID SIMULATION 

A concept of hybrization between the MD simulation and 
the ACVT simulation is that the ACVT simulation works when 
a kinetic energy of a projectile is higher than the threshold 
energy £^th = 200 eV, while the MD simulation works when 
the kinetic energy becomes lower than the threshold energy 
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FIG. 1: Schematic diagram of the ACVT-MD hybrid simulation of 
hydrogen injection onto a graphite. 

To explain the hybrid simulation more concretely, let's con- 
sider a hydrogen atom injection onto a graphite (Fig. [T]). An 
injection kinetic energy of the hydrogen atom is set to 1 keV, 
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FIG. 2: Data flow diagram between the processes of the ACVT-MD 
hybrid simulation code. 



which is higher than E'th . While the kinetic energy of the hy- 
drogen atom is higher than £^th, trajectories of the hydrogen 
atom and the surrounding carbon atoms are calculated by the 
ACVT simualtion. The hydrogen atom dissipates its kinetic 
energy by interacting with carbon atoms in diamond and then 
the kinetic energy becomes lower than £^th. At that moment, 
the MD simulation starts to calculate the motions of the hy- 
drogen and the carbon instead of the ACVT. Some of carbon 
atoms which collide with the incident hydrogen are kicked out 
on the way of the ACVT simulation. The motions of all kicked 
carbons in the cascade process are solved by the same proce- 
dure as the incident hydrogen. 



In order to realize the above hybrid simulation, we use Mul- 
tiple Program Multiple Data (MPMD)^ which was developed 
to execute different programs parallely in one simulation with 
synchronizing their data. Using the MPMD, we control four 
programes, i.e., (i) master program, (ii) ACVT program, (iii) 
MD program and (iv) splitting program as shown in Fig. O 
The master process is installed between the ACVT and MD 
processes for data synchronization. After an initialization of 
each process (Fig.[2l) and (2)), the ACVT process (3) is exe- 
cuted until the kinetic energy of the incident hydrogen and all 
kicked atoms become lower than £^th- When the kinetic ener- 
gies of all atoms becomes lower than £^th, the ACVT process 
sends, to the master process, a list of the position and velocity 
of the moving atoms as a result of the ACVT simulation. The 
master process (4) informs the splitter process of the comple- 
tion of the ACVT process with the list of moved atoms. The 
splitter process updates its material data with the list of moved 
atoms. (5) Then, the MD simulation box are picked up from 
the entire target material and they are sent to the MD pro- 
cesses. The MD simulation is started in each MD processes. 
(6) We will explain the algorithm of the ACVT simulation, 
the MD simulation and the splitting program in the following 
subsections. 



A. ACVT Simulation 

A binary collision between a projectile and a target atom is 
necessary in the ACVT simulation*^ . During the binary colli- 
sion, a total energy, a total momentum and a total angular mo- 
mentum are conserved. Due to these conservations, the final 
position and velocity of the projectile and the target atom are 
obtained analytically on a potential V(t), where r is the dis- 
tance between the projectile and the target atom. In the ACVT 
code similarly to the AC AT code, the Moliere approximation 
is adopted and then Vir) is given by the Thomas-Fermi po- 
tential as follows: 
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where a is a screening length, Zi and Z2 are atomic num- 
bers of the projectile and the target atom, respectively. Unlike 
the MD simulation, the ACVT simulation directry gives us the 
asymptotic trajectory of the projectile collided with the target. 
In the ACVT simulation, the projectile hits the target atoms 
and then it brings about cascade shower in the large simula- 
tion box. 



B. Molecular Dynamics Simulation 

The equation of motion in the MD simlation is written by 
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where ri , pi and rrii are the position, momentum and the mass 
of the i-th atom. The function U is total potential energy. In 
our MD simulation, the modified REBO potential is used for 
the potential U which is defined as: 
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where rij is the distance between the i-th and the j-th atoms. 
The functions V^^^ and V^^.^ represent repulsion and attrac- 
tion, respectively. The function bij generates a multi-body 
force. The bond angle Oj-j^ is the angle between the vector 
from the i-th atom to the j-th atom and the vector from the 
i-th atom to the k-th atom. The dihedral angle O^^i is the 
angle between the plane passing through the k-th, i-th, and 
j-th atoms and the plane passing through the i-th, j-th, and 
l-th atoms. Moreover, the parameters {r}, {0^}, and {6'^^} 
denote all sets of rij , Of-j^ , and 6^^i , respectively (for details 



of the modified Brenner REBO potentialm, see Refs*^ andiS). 
The second-order symplectic integration^^ is used to execute 
the time integration of the equation of motion where the time 
step is set to 5 x 10~^^ s. 

Unlike the ACAT simulation, the MD simulation cannot 
treat a large scale material. To calculate efficiently, local boxs 
14 A on a side are picked up, and the simulation for the atoms 
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FIG. 3: Structure of a MD simulation box. 
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in each box is executed in parallel. To obtain precise forces 
acting on all atoms in the MD simulation boxs 14 A on a side, 
the bond energy Uij in Eq. ^ is calculated in a larger cubic 
box 26 A on a side. Furthermore, the function bij depends 
on more outer covalent bonds. Consequently, the MD pro- 
cess needs the positions of atoms in a cubic box 34 A on a 
side. Namely, in the MD process, the equations of motions 
are solved about atoms in the cubic box 14 A on a side, and 
forces acting on the atoms are calculated by aggregating sur- 
rounding atoms in a cubic box 34 A on a side as shown in Fig. 

m 



C. Detamination of MD Simlation Box 

By the MPMD method, one CPU is allocated into the MD 
process treating one cubic box 34A on a side. However, there 
is the special case that two or more cubic boxs overlap be- 
cause the moving particles are close. These overlapping boxs 
are regarded as a combined MD simulation box for one CPU 
process. The determination for the MD simulation box is car- 
ried out by the splitting program (See Fig.|4l) as follows. The 
simulation box in the ACVT code is divided into cubic cells 
2 A on a side, where 2 A is the cutoff length of covalent bond 
in the REBO potential. 17x17x17 cells, which are corre- 
sponding to the cubic box 34 A on side, are marked up. Only 
a center cell of the marked cells is a member of the MD sim- 
ulation box, initially. The center cell is set to the first target 
cell and the member of the MD simulation box is selected by 
the following recursive routine: if adjoining cells of the target 
cell has been marked up, they are added into the member of 
the MD simulation box and become next target cells, while if 
they have not been marked up, the recursive routine is stopped. 



III. NUMERICAL DEMONSTRATION 

We demonstrate an ACVT-MD hybrid simulation for a hy- 
drogen atom injection onto a graphite, whose size is 321 A 



FIG. 4: Example of the determination procedure of the MD simula- 
tion box. The gray and white cells indicate the marked cells as MD 
region, respectively. The black cells indicate the counted cells by the 
recursion routine. (l)Dividing the target material into cell. (2)Mark- 
ing the cells corresponding to three target atoms. (3),(4)Counting the 
MD simulation box A by recursion routine. (5)-(9)Counting the MD 
simulation B which includes two target atoms. 



TABLE I: The average of the depth d± and the horizontal distance d\\ 
of 900 injections. 9% trajectories which go out of the target material 
are excepted. 



[A] 



[A] 



ACVT-MD 127.3 ± 31.5 59.1 ±34.3 
ACVT 127.2 ± 37.5 82.13 ± 45.8 



long, 347 A wide and 335 A depth which consists of 100 
graphene. Periodic boundary conditions are imposed on the 
horizontal direction. The initial temperature of the graphite 
is set to zero Kelvin. The hydrogen atom is injected normaly 
onto the surface at 1 keV. 

900 simulations with the same initial material and ran- 
domly changed injection position have been performed. As 
a bench-mark, stand alone ACVT simulations have also been 
performed. Table II shows the average and the standard devi- 
ation of the depth d± and the horizontal distance from the 
surface to the final positions in the result of both the ACVT- 
MD hybrid simulation and the ACVT simulation. Although 
the averages of depth in the both cases are almost equal, the 
average and the standard deviation of horizontal distance by 
the ACVT-MD hybrid simulation are smaller than that of the 
ACVT simulation. Figure Oa) and (b) show 30 samples cho- 
sen randomly in 900 trajectories calculated by the ACVT and 
the ACVT-MD hybrid code, respectively. Figure Ob) shows 
the trajectories in the same initial condition calculated by the 
ACVT-MD hybrid simulation. The gray and blue lines in 
Fig. [5tb) show the trajectories of which indicate the ACVT 
contribution and the MD contribution, respectively. Compar- 
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FIG. 5: 30 sample trajectories of incident hydrogen that move in 
graphite calculated by the (a) stand alone ACVT code and (b)ACVT- 
MD hybrid simulation, respectively. (c)Detailed image of the MD 
simulation part for trajectory 'A' in (b). 



ing these two figures, it is found that the ACVT-MD hybrid 
simulation gives shorter trajectories than the ACVT simula- 
tion does. As shown in the trajectory 'A' in Fig. [3a), the 
incident hydrogen atom moved the long distance along the in- 
terlayer region of the graphite, because the ACVT code does 



not treat the binding from surrounding atoms. The MD simu- 
lation includes it by the multi-body force. Therefore, the atom 
is not possible to move a long distance for low kinetic energy 
in the ACVT-MD hybrid simulation as shown in Fig.Oc). In 
our previous research^, it was shown that there is an adsorp- 
tion site 1.1 A above each of the carbon on a graphene sheet. 
The hydrogen calculated by the MD dissipates its kinetic en- 
ergy quicker than the ACVT case and finally be trapped at the 
adsorption site. 

IV. SUMMARY 

The ACVT-MD hybrid simulation code has been developed 
to treat large materials, long time and wide range kinetic en- 
ergy in PSI. The motion of hydrogen atom is solved by BCA 
method using the ACVT code while its kinetic energy is higher 
than 200 eV. When the kinetic energy becomes lower than 
200 eV, the MD code is invoked around atoms moved by the 
ACVT simulation. 

The ACVT-MD hybrid simulation was performed and com- 
pared with stand alone ACVT simulations of 900 hydrogen 
atom injections onto a graphite. Their results suggested that 
the ACVT-MD hybrid simulation tends to derives a shorter tra- 
jectory than the ACVT simulation does. This fact is caused by 
the binding from the surrounding atoms, which is ignored in 
the ACVT simulation. 
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